Electron correlation in the Si(100) surface 
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Motivated by the controversy between quantum chemists and solid-state physicists, and by recent 
experimental results, spin-polarized density-functional (DFT) calculations are used to probe electron 
correlation in the Si(100) reconstructed surface. The ground state displays antiferromagnetic spin 
polarization for low dimer inclinations indicating, not magnetic order, but the importance of Mott- 
like correlations among dangling bonds. The lowest energy corresponds to a higher dimer inclination 
with no spin. DFT energies, however, should be taken with caution here. Our results together with 
quantum-chemical findings suggest dimers with highly correlated electrons that tend to buckle due 
to interactions with other dimers. 
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The freshly cut (100) surface of silicon has two dangling 
bonds per Si atom. This is a highly unfavorable con- 
figuration, and the surface reconstructs forming rows of 
dimers, thereby reducing the number of dangling bonds 
by a half. This chemical description is clear and well es- 
tablished (jjj, but the physics of the remaining two dan- 
gling bonds per dimer is not. The main reason is that 
the CT-like chemical bonds in silicon are well described 
by many electronic-structure methods, whereas the dan- 
gling bonds are much more subtle. The early argu- 
ment for asymmetric dimers based on relating symmetric 
dimers to metallic character Q was found to be incorrect, 
since it is possible to obtain a symmetric, non-metallic 
state within a Hubbard like framework |3] ||. However, 
experimental evidence accumulated afterwards strongly 
suggesting the asymmetry of the dimers. In addition, 
solid-state first-principles calculations based on the local- 
density approximation (LDA) to density-functional the- 
ory (DFT) converged to results in agreement with ex- 
periments |(|. The dimer asymmetry was rationalized 
in chemical terms, the driving force staying within each 
dimer: a rehybridization of the orbitals with a charge 
transfer from the lower to the upper atom, leaving the 
lower with an sp 2 -like hybridization that would open the 
one-particle gap and thus lower the total energy |l],f7j . 

The issue about the symmetry /asymmetry of the 
dimers would seem to be settled and of no further rele- 
vance now, except for two facts: (i) Kondo et at. || have 
recently performed low-temperature (20 K) scanning tun- 
neling microscopy, and they found that at these low tem- 
peratures the "symmetric dimers dominate the surface," 
with the asymmetric dimers being dominant at higher 
temperatures (110 K). Based on these results they claim 
that the ground state for the Si(100) system consists of 
symmetric dimers, and that the theory should be revised 
in the line of correlated electrons. It is worth mentioning 
that a similar result was also obtained by another ex- 
perimental group using the same technique M, but they 



propose that the reason for the observed symmetric ap- 
pearance of the dimers is an anomalous flipping motion of 
the buckled dimers; and (ii) different groups in quantum 
chemistry [|l0|— |l2|] have shown that a proper treatment of 
electron correlations produces symmetric dimers with a 
highly correlated ground state [HJi). 

From the point of view of the theory, the controversy 
can be explained by noting that: (i) solid-state DFT cal- 
culations are performed for realistic infinite-surface (slab) 
geometries, but cannot describe the electronic correla- 
tions in a controlled way; and (ii) quantum-chemistry 
calculations can describe electron correlations properly, 
but have to work on isolated clusters containing up to 
a few dimers. It is clear that a theory with appropriate 
geometry and correlation description is needed. A the- 
ory that (we conjecture) would show that if the dimers 
are asymmetric it is because of reasons extrinsic to the 
dimer, i.e., interactions among dimers. The dimers would 
display important Mott-like (antiferromagnetic) correla- 
tions, and there would be no driving force towards asym- 
metry within the dimer. 

Spin polarization is the simplest way to assess the 
importance of Mott-like correlations. A spin-polarized 
wave-function tends to the right dissociation limit, ex- 
cept for the fact that the spin symmetry is artificially 
broken. A spin-polarized DFT study of the surface is 
presented below. The polarization itself is to be inter- 
preted as an indication of the need of a better treatment 
of dynamical correlations, and, therefore, the numbers 
are to be taken with caution. The magnetic order is arti- 
ficially imposed and of no relevance to the point, except 
as a way of describing the short range correlations. 

The DFT calculations were made using the numerical 
atomic-orbital (NAO) method (13|hT5| in the Siesta code 
implementation The generalized gradient approx- 

imation (GGA) of Perdew, Burke, and Ernzerhof Jlfj] 
to Kohn-Sham jlTj theory was the functional for ex- 
change and correlation, except for the LDA JTJ] tests 
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shown below. Core electrons were replaced by norm- 
conserving pseudopotentials fjlJ in their fully non-local 
formulation |20|. Valence electrons were described using 
a double-^ polarised basis set of numerical atomic orbitals 
with a range defined by an energy shift of 200 meV, and a 
variationally optimized splitnorm of 0.3 [ fl5| . A uniform 
real-space mesh with a plane-wave cutoff of 80 Ry was 
used for numerical integration. 

Calculations were performed for the 2x1, p(2x2) and 
c(4x2) reconstructions of Si(100) using a repeated slab 
geometry with nine layers of silicon and a layer of hydro- 
gen atoms saturating one of the sides of the slab. The 
H-Si bonds were relaxed in a bulk-like saturated slab with 
H on both surfaces, using the bulk lattice parameter ob- 
tained within the theoretical framework described above. 
Then, one of the surfaces was kept saturated and fixed, 
including the H layer and the closest Si layer. The re- 
maining eight layers were allowed to relax freely in all 
the calculations described below. Relaxations were con- 
sidered to be finished when the maximum residual force 
was below 0.02 eV/A. The forces on the fixed layers al- 
ways stayed <0.06 eV/A. Integrations over the first Bril- 
louin zone were approximated by Monkhorst-Pack ]2l| ] 
sets of fc-points for a length cutoff of 15 A p2| . which 
corresponds to a 4x8x1 set in the 2x1 geometry. 

To test the methodology we have reproduced earlier re- 
sults obtained with LDA. For bulk silicon a lattice param- 
eter of 5.40 A is obtained, which compares very well with 
the experimental value of 5.43 A [^3| and with 5.37-5.38 A 
of LDA plane- wave (PW) calculations The bulk 

modulus is 95 GPa employing the Murnaghan equation 
of state p^J and 96 GPa using a quartic fit, to be com- 
pared with 99 GPa of experiment j2^| and 96-98 GPa of 
PW [ p4p5| ]. Table I shows the relevant magnitudes ob- 
tained for the 2x1, p(2x2) and c(4x2) reconstructions 
as compared with the figures obtained from LDA calcu- 
lations of other groups. Spin polarized results for LDA 
were never reported, probably because the calculations 
converged to non-polarized. We obtain non-polarized re- 
sults within LDA. 

The spin-polarized GGA results of the present study 
are summarized in Fig. 1 for the p(2x2) geometry, where 
the relevant magnitudes are plotted versus dimer buck- 
ling. The buckling angle was imposed as a constraint in 
each calculation, with the rest of the system's degrees of 
freedom being free to relax. As can be observed from 
Fig. 1, the most stable configuration remains unpolar- 
ized, as in the LDA calculation, with asymmetric dimers 
of 18.5° of inclination. At small buckling angles (<6.1°) 
the system has a polarized, antiferromagnetic (AF) or- 
dered solution with a lower energy than the non-polarized 
solution (0.07 eV/dimer decrease at 0°). E(a) seems to 
show a shallow minimum around the symmetric config- 
uration. The dimer bond length increases with the spin 
polarization. The shorter bond length of 2.32 A of the 
unpolarized solution can be understood in terms of the 



TABLE I. LDA results. Dimer bond length (d), buckling 
angle (a), energy difference with respect to symmetric 2x1 
relaxed configuration (AE), displacement of subsurface atoms 
in the direction parallel to the dimer rows (Ax), and surface 
energy (E s ) M. 
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(artificial) reinforcing of the bond (double-bond) in the 
dimer, essentially absent in the polarized solution, for 
which a typical Si bond length of 2.35 A is obtained. 
The subsurface lateral relaxation correlates with the in- 
clination, but also with the dimer bond length. 

A magnetic moment of 0.47 \ib is obtained at 0° in- 
clination for this AF configuration. For a 2x1 unit cell 
there are two additional spin-polarized solutions: a fer- 
romagnetic (F) spin configuration, and a mixed one, AF 
within the dimer and F among dimers in a row. Their en- 
ergies are higher than the AF and are shown in Fig. 1(a). 
The magnetic moments are 0.51 hb and 0.35 (1b for the 
mixed and F configurations, respectively. Interestingly, 
the dimer bond lengths are 2.31 A and 2.33 A for the 
mixed and F configurations, respectively, shorter than in 
the AF case. 

These results confirm the importance of Mott-like elec- 
tronic correlations in this surface. The results do not 
tell, however, about the intrinsic or extrinsic origin of the 
buckling that appears. Intcrdimer interactions that can 
make a difference are: (i) electrostatic: each dimer has 
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FIG. 1. GGA results versus dimer inclination, a, with and 
without spin polarization. The curves are interpolations. 

(a) Total energy (per dimer) with respect to the most sta- 
ble p(2x2) configuration. The full circles are for the ground 
state, which is spin polarized (AF) at small angles; empty 
circles show the unpolarized solution. The square is for the 
mixed AF-F configuration described in the text, and the dia- 
mond for the ferromagnetic configuration. 

(b) Magnetic moment on each dimer atom in the AF config- 
uration. 

(c) Intradimer charge transfer (Qhigh ~ Qiow), obtained from 
Mulliken population analysis. 

(d) Dimer bond length. 

(e) Displacement of subsurface atoms in the direction parallel 
to dimer rows. 



an electric dipole due to its charge transfer p| ]34|]35H an d 
a polarizability that can enhance or depress the dipole 



when interacting with its neighbors; (ii) surface stress: 
the relative disposition of dimers can relax surface stress 
as pointed out by Garcia and Northrup [Q; (in) elec- 
tronic hopping: the electron derealization among dan- 
gling bonds of different dimers can depress electronic cor- 
relations and change the energy. 

The first point can be discarded by noting that it has 
opposite effect on the energy for the asymmetric 2x1 and 
p(2x2) surfaces. The difference between these two struc- 
tures is in the relative disposition of neighboring dimer 
inclinations: in the 2x1 all the dimers buckle in the same 
direction, in the p(2x2) neighboring dimers in a row have 
opposite inclinations. The electric dipoles are mostly par- 
allel and repelling in the 2x1 order, and antiparallel and 
attracting in the p(2x2). The effect should thus appear 
in the energy difference between the asymmetric 2x1 and 
thep(2x2) rather than in the difference between symmet- 
ric and asymmetric. This argument has been ratified by 
calculating the electrostatic interaction energy of a two- 
dimensional set of dipoles with the appropriate disposi- 
tion and the dipolar moments extracted from the DFT 
calculations. The results confirm the order of magnitude 
and sign used in the argumentation: the 2x1 disposition 
has positive energy, whereas the p(2x2) has it negative, 
both having an order of magnitude of 0.1 eV/dimer for 
the maximum p(2x2) charge transfer of 0.36 electrons. 

It is tempting to use a similar argument for the effect 
of the surface stress: it seems reasonable to assume that 
by buckling the dimers the surface stress would behave 
in opposite manner for the 2x1 configuration than for 
the p(2x2) one. However, the calculations do not show 
that trend: the compressive stress in the direction of the 
dimer rows is highest for the symmetric-dimer surface, 
and lowest for the p(2x2) structure. Surface stress could 
thus be playing a role in the destabilization of the sym- 
metric dimers. 

The effect of the electron derealization among dimers 
cannot be assessed as easily. The subsistence of im- 
portant correlation effects can be inferred from the fact 
that clusters with more than one dimer, when described 
within quantum-chemical methods that describe well the 
dynamical correlations, tend to stay symmetric |0,B6[. 
However, to our knowledge, no high-level quantum- 
chemical calculation comparing the buckled against the 
non-buckled configurations, has been performed for clus- 
ters with more than two surface dimers. As has been 
shown recently by Penev et al. ]37j this may be an im- 
portant point, as discussed below. Penev et al. |57j help 
in discerning the issue of electron derealization by com- 
paring cluster calculations with slab calculations within 
the same level of theory otherwise. The asymmetry of the 
dimers is found to be very much favored in the extended 
system, the clusters with one or two surface dimers A 
showing a nearly flat behavior of energy versus buckling 
angle, E(a), close to a = 0°. Part of the difference, 
however, could be due to the treatment of the relaxation 
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constraints in the clusters, as pointed out by Yang and 
Kang fjSf . An approximation to the correlated and ex- 
tended situation could be obtained by 

E{a) » E%£(a) + {E°[ b T (a) - E^{a% (1) 

where E^°^(a) is for the correlated cluster calculations, 
and E^ b T (a) and E^^(a) are for the slab and clus- 
ter DFT calculations, respectively. From the calcula- 
tions of Penev et al. |37]] one can estimate that the term 
[E^ b T (a) — E®^(a)] would give a maximum gain of 
the order of 0.2 eV/dimer towards buckling due to ex- 
trinsic dimer effects (comparing the slab with the SigHi2 
cluster calculation). To determine if the ground state has 
a buckled or non-buckled configuration one has to know 
^cius( a ) as a function of a for a high- level correlated 
calculation. A more rigorous computation for the corre- 
lated extended system requires computation techniques 
not available at present. Given all that was discussed 
above, a DFT calculation, spin-polarized or not, should 
not be fully trusted in this respect. 

Another possibility, which can also help to gain fur- 
ther insight about the correlation effects in this problem, 
would be the use of a model Hamiltonian for a single 
row of dimers, as the stripe or ladders that are being 
treated now for High T c superconductors, as pointed out 
by Kondo et al. || . It would consist of a Hubbard-like 
model on the ladder, coupled to a classical dynamical 
variable describing the buckling. Parameters for this 
model can be obtained from suitable calculations with ac- 
cessible techniques, and the model can be solved in differ- 
ent approximations to a high degree of accuracy. Results 
for it will be presented elsewhere (39|. It is interesting 
to note that such a Hamiltonian is remarkably close to 
Hamiltonians used for completely different problems in 
many-body physics. There are thus cross-breeding pos- 
sibilities, and it would be quite interesting if knowledge 
on the Si(100) surface could be used to gain insights for 
the correlated motion of carriers in doped perovskites. 
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